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Abstract 

A general two-dimensional Euler zonal method 
has been developed for computing flows about com- 
plex airfoil geometries such as multielement and 
Iced airfoils. The method utilizes a coifiposite 
structured and unstructured grid generated using 
conformal mapping and Delaunay triangulatlon, 
respectively. The finite-volume Euler method Is 
then modified to couple solutions In the zones 
with structured and unstructured grids. Solutions 
about an Iced airfoil and a multielement airfoil 
are given as examples of applications of the 
scheme . 

1 .0 Introduction 

The aerodynamic analysis of complex airfoil 
sections continues to receive much attention In 
experimental, theoretical and computational stud- 
ies. Such analyses of airfoil flows Include Iced 
airfoils, multielement airfoils, and advanced air- 
foil concepts with divergent trailing edges. Con- 
ventional computational methods based on a single 
zone mesh, developed for simple airfoil configu- 
rations, are generally not suitable for these 
geometries since a single structured grid of suf- 
ficient quality cannot be generated. Zonal or 

unstructured-mesh methods that can provide ade- 
quate mesh resolution near high pressure gradient 
regions In each zone are needed to handle the 
complexity of the flowflelds. An attractive zonal 
approach 1s to use structured meshes In most parts 
of the domain and unstructural meshes In enclosed 
regions next to the portions of airfoil that are 
difficult to model with structured grids. The 
objective of this paper Is to develop a general 
two-dimensional Euler method based on this 

approach. 

Although many methods have been developed for 
the generation of structured grids around simple 
geometries, few can be readily extended to comp- 
licated two- and three-dimensional shapes such as 
multielement or Iced airfoils or multicomponent 
aircraft configurations. Also, In many applica- 
tions of structured grids, the quality of the 
generated mesh Is not uniform throughout the comp- 
utational domain, resulting In poor resolution In 
specific regions of Importance. For this reason, 
triangular (tetrahedral) meshes have proven very 
attractive, since relatively complex geometries 
can be meshed efficiently, and an almost arbitrary 
degree of mesh adaption and refinement can be 
achieved by the addition of control points. Fur- 
thermore, the development of flow algorithms that 
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do not depend on the Inherent structure of the 
grid points have eliminated the restriction of 
grid structure and made triangles and tetrahedrons 
suitable shapes with which to resolve complicated 
geometric regions. 

The debate over relative advantages of struc- 
tured and unstructured grid methods Is an ongoing 
matter. While unstructured grid methods have the 
clear advantage In that they can treat complicated 
problems such as a complete aircraft configuration 
under static or dynamic deformation,^ they have 
been regarded as less efficient and less accurate 
than their structured counterparts .2.3 Lg^k of 
robust acceleration techniques, such as multigrid 
schemes for three-dimensional problems. Is also a 
clear disadvantage. Furthermore, direct and 
Implicit solvers for unstructured grids are devel- 
oped, but In general are not as efficient and 
effective as their structured counterparts. 
Structured flow solvers have many other highly 
desirable features Including efficient grid gen- 
eration techniques and smaller computer time and 
memory requirements. 

An efficient way of analyzing a complex geom- 
etry of several components Is through a zonal 
approach, using composite structured and unstruc- 
tured grids. This approach requires considerably 
less memory than using an entire unstructured mesh 
capable of handling the same geometry. In view of 
the different desirable features of structured and 
unstructured meshes, the present zonal approach 
takes advantage of structured meshes In appropri- 
ate regions while using the versatility of the 
unstructured grids In others. 

In this paper, a general two-dimensional zonal 
boundary Interactive scheme utilizing combinations 
of structured and unstructured mesh types Is pre- 
sented, Solutions for flow about multielement and 
Iced airfoils are given as an example of the 
application of the scheme. 

2.0 Grid Generation 

A brief grid generation process will be des- 
cribed In this section. A base structured grid 
which encompasses the entire flowfleld Is first 
generated using existing methods. Regions of 
undesirable grid quality are then Identified and 
categorized as subsequent unstructured-mesh zones. 
Triangular grids can be generated In these regions 
using existing mesh points and additional points 
as necessary . An unstructured grid generation 
method based on Delaunay triangulatlon, developed 
earlier, is applied for this purpose. An 

essential requirement for obtaining satisfactory 
meshes when using the Delaunay triangulator Is 
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appropriate placement of Interior points, since 
this scheme gives no guidance on where to place 
the mesh points. Inappropriate placement of mesh 
points results In poor quality meshes, even though 
they satisfy the Delaunay criteria. Therefore, 
certain criteria for Introducing Interior mesh 
points have to be established before applying the 
Delaunay trlangulatlon. 

The Delaunay trlangulatlon method Is applied In 
three steps to generate triangular meshes In the 
unstructured zone enclosed by the surrounding 
structured meshes and solid surface boundary. The 
first step Is to triangulate the points along the 
body and zonal boundaries In order to obtain an 
Initial trlangulatlon based solely on existing 
boundary points. 

The second step deals with the placement of 
Interior points. It Is Important to place a suf- 
ficiently dense mesh of points In high gradient 
regions such as corner regions, leading- and 
tralling-edge regions, etc. A combination of 
C-mesh points around the leading edge and 
Cartesian mesh points In other parts pf the 
domain can be chosen as the Interior points. 
Based on Bowyer's algorithm, a series of new 
points are added one by one to the existing trl- 
angulatlon. removing triangles close to the point 
being Inserted, and reconnecting the new point to 
the existing nodes In such a way as to form new 
triangles which satisfy the Delaunay criteria. 
This procedure Is repeated for all new points 
Introduced In the domain. 


The final phase Is a grid-smoothing procedure. 
This Is necessary because points that are Intro- 
duced sometimes fall too close to each other, 
resulting In skewed triangles, such as triangles 
with large aspect ratios. An Iterative Laplaclan 
smoothing scheme Is applied to Improve the overall 
quality of the unstructured mesh. 

3.0 Finite Volume Scheme and T1me-Stepp1n_g 

Finite volume Euler methods of Jameson® and 
HavrlpHs and Jameson^, using fourth-order Runge- 
Kutta time-stepping, are adopted to Interactively 
solve for the flows In different zones. These 
methods are briefly reviewed here to point out 
some modifications that are made 1n the present 
work . 

The Euler equation for two-dimensional 1nv1sc1d 
flow In integral form for a region n with a 
boundary 3J5 Is given as 


fr Jj wdxdy ♦ J (fdy - gdx) * 0 (1) 


where x and y are Cartesian coordinates and 
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Here p, u, v, p, E and H are density, velocity 
components, pressure, total energy, and total 
enthalpy, respectively. For a perfect gas we have 
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Thus, we only need to solve for four variables: p, 
pu, pv and pE. Equation (1) Is discretized over 
Individual control volumes (triangles) In a cell- 
centered approximation In which the flow variables 
are stored at the center of each cell. The above 
discretization procedure with the addition of 
artificial dissipation terms results In a set of 
ordinary differential equations. 

dw. 

Si t [Q(w^) . D(w^)] . 0 (3) 

where Is the area of the cell. Q Is the 
spatial approximation of fluxes given by the sec- 
ond part of Eq. (1), and D Is an appropriately 
constructed dissipation operator. 


The fourth-order Runge-Kutta scheme Is used to 
advance the solution In time from time step n to 
time step n+1. With the nonlinear operator P 
defined as 

P(W) . 1 [0(W) - D(w)] (4) 

we have 
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The artificial dissipation operator Is calculated 
In the first and third step only to save computing 
costs. 


In the region with structured grids implicit 
residual averaging Is used to Increase the base 
of the CFL number. It Is shown'® that Implicit 
smoothing increases the stability to Courant num- 
bers much greater than the Courant number limit 
of the explicit scheme. The calculations pre- 
sented In this work are obtained with a CFL number 
of 7. In addition, a variable time step, based 
on the maximum stability limit set by the local 
Courant number, and enthalpy damping are used to 
accelerate the convergence of the solution. Dis- 
sipation terms are formed as a blend of second- 
and fourth-order terms, coefficients of which are 
adapted to the flow. The resulting scheme 1s 
second-order accurate 1n the regions of smooth 
flow and first-order accurate near the shocks. 
To ensure that these dissipative terms are sig- 
nificant only 1n the vicinity of a shock, the 
second-order dissipation terms are scaled by the 
local Laplaclan In the pressure. 

Similar to the structured flow solver, In the 
unstructured flow solver fourth-order Runge-Kufta 
time stepping 1s also used to advance the solution 
In time. The support of the time-stepping scheme 
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here Is, however. Increased by an explicitly 
residual averaging scheme. If the residuals at 
cell 1 are 
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The second-order dissipation term Is similarly 
constructed by replacing and In Eq. ( 6 ) 

with and W|^. respectively, so that we have 
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where are residuals at three forming points 
of the triangular cell. The residuals at nodal 
points are obtained as the average of residuals 
at all cells having that point In common, c Is 
a constant which Is chosen as 0.6. With this 
smoothing scheme, the CFL number for the unstruc- 
tured flow solver could also be Increased to 
about 7, 4 

It Is found that the CFL number In the two 
zones must be similar to facilitate the converg- 
ence of the solution In the entire domain. 

4.0 Dissipative Terms 

The structured flow solver uses a blend of 
second- and fourth-order dissipation terms to 
prevent odd and even decoupling of the solution. 
The artificial dissipation for the unstructured 
flow solver Is similarly constructed as a blend 
of undivided Laplaclan and biharmonic operators. 
Generally, 1n the absence of shocks In subsonic 
flows, only biharmonic dissipation Is required. 
However, It Is found that In regions of large 
pressure gradients such as at the leading edge of 
the flap, second-order dissipation terms are still 
needed even In the subsonic flow regime. 

To obtain the fourth-order dissipation term In 
the triangular mesh zone, an undivided four-point 
Laplaclan operator 1s first defined as 

2 3 

w = I w - 3w 
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where w represents the flow variables p, pu, pv, 
and pH. The dissipation flux across a cell face, 
Ik, delimiting cells. 1, and Us neighbors, k, are 
then calculated as 

‘'ik • 

where 
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Here and ayj^ are coordinate Increments of the 
edge and uj^, vj^, Cj^ are velocity components and 
the speed of sound along the edge, respectively, 
and are taken as the average of the values at 
cells 1 and k. The A|^ term, which Is propor- 
tional to the size of the cell face, k, and repre- 
sents the maximum eigenvalue of the Euler equa- 
tion In the direction normal to the face, scales 
appropriately with the time derivative In Eq. (3). 
In order to more accurately scale the cells which 
have higher aspect ratios, 1s not Integrated 
around the boundary of ,y\e control volume, as Is 
commonly done. The c 1s a constant defined 
later. 


where should be of order one near a shock and 

and of order (ax)^ 1 n regions of smooth flow to 
preserve the secopdrorder accuracy of the scheme 
To ensure this. Is then scaled proportional 
to an undivided Laplaclan In the pressure 
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It Is found, however, that a better scaling factor 
for the edges that are on the solid boundary Is 
the pressure gradient In the surface-wise direc- 
tion. Thus, the pressure gradient term for cells 
on the solid boundary Is defined here as 


•*1 = < 


PUI - * Pj-1 

Pui " 2Pi ♦ Pl-1 


) 


where 1 , 1 - 1 , and 1+1 are the adjacent cells to 
cell 1, which are on the surface boundary. For 
these cells, cf|| Is then taken as 

( 2 ) 


It Is known that the fourth-order dissipation 
term may produce overshoots In the vicinity of a 
shock; therefore, they are turned off by defining 
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Here c 2 and C 4 are empirically determined 
constants . 


The final 
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form of the dissipative flux Is then 
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This dissipation flux, d^l^, Is added or sub- 
tracted from Its adjacent cells consistent with 
the direction of the normal to each edge. This 
ensures that dissipation 1 s conserved throughout 
the region. 


The treatment of the boundary conditions Is 
known to affect the convergence and accuracy of 
the solution. The wall pressure can be extrapo- 
lated from the pressure at the center of the 
boundary cell using the method given In Ref. 8 . 
However, In the present calculations It Is found! 
In both the structured and unstructured solvers! 
that the normal pressure gradient, ap/an. Is 
negligible, and the pressure of the boundary could 
be assumed equal to the pressure at the center of 
the boundary cel 1 . 


It Is Important to note that when the solid 
boundary condition Is Imposed In the unstructured 
flow solver, however, Is that setting the normal 
components of the fluxes to zero and only account- 
ing for the pressure terms In the momentum equa- 
tion does not necessarily satisfy the flow tan- 
gency condition on the boundary. A stronger 
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formulation is needed to ensure this requirement. 
In the present cell-centered scheme, a stronger 
form of the boundary condition Is Imposed In order 
to compute the artificial dissipation terms more 
exactly. Flow variables #>, pu, pv, H In Imaginary 
cells Inside the solid surface are extrapolated 
based on the assumptions of no normal flux and 
equal tangential fluxes between cells outside and 
Inside the boundary. These values are then used 
In the calculation of dissipation terms associated 
with the edges that are on the boundary. This can 
be shown to be equivalent to explicitly setting 
velocities to be tangent to the wall, as suggested 
in Ref. 11. 


5,0 Zonal Interference Scheme 

Interaction between neighboring zones could be 
greatly simplified by ensuring that neighboring 
structured and triangular grids have comf>lete, 
common, edges. This could easily be done alon^ 
the zonal boundaries by choosing the structured 
grid points as forming points of the triangular 
grids. Mass, momentum and energy fluxes through 
the zonal boundaries are conserved by using the 
values from neighboring zones when calculating 
fluxes. This requires maintaining Information on 
the grid Interaction, Including Indexes of struc- 
tured grids next to the boundary triangular 
meshes, and vice versa. 

Careful attention must be given to calculation 
of dissipation terms for the boundary cells In 
order to ensure conservation of dissipation terms 
throughout the flow. Dissipation terms for the 
triangular, boundary cells are calculated by 
extracting Information from both the structured 
and unstructured . zones and vice versa. Inappro- 
priate treatment of dissipation terms can result 
in solution Inaccuracy along the zonal boundaries 
or can produce considerable decoupling of the 
solution In boundary regions of both zones and 
contaminate the solution In the entire region. It 
is also found that comparable CFL numbers and 
degrees of smoothing In different zones, are 
essential to Improve the overall convergence of 
the solution in the entire region. 

6,0 Results 

The multielement airfoil of Fig. 1 with over- 
hang of 514, gap-to-chord ratio of 1014, and flap 
angle of rotation of 15® is considered as the 
first computational example. The Isobar solution 
at a Mach number of 0.2 1s shown 1n F1g. 2, which 
indicates a smooth solution across the zonal 
boundaries. This solution Is obtained by using 
256 X 32 0-meshes, part of which are replaced by 
976 triangular meshes. A CFL of 7, along with a 
variable time step based on a maximum limit set 
by local Courant number, Is used In both zones. 

Figure 3 shows the details of the solution In 
the unstructured zone and indicates that the 
strong form of the boundary condition discussed 
earlier In Section 4.0 adequately satisfies the 
flow tangency requirements on the solid boundary. 
The solution also accurately pred1cts_the location 
of the stagnation point on the leading edge of the 
flap. Furthermore, even though the equations are 
Inviscid, the addition of dissipation terms Intro- 
duces viscous flow like vorticity In the flap 
well. 


Calculated surface pressures for both the main 
airfoil and the flap are compared with experi- 
mental data of Ref. 12 In Fig. 4(a,b). The exper- 
imental results are obtained at a M > 0.195 and 
Reynolds number of 500,000 based on unextended 
chord length. The comparison 1s obviously very 
poor In the flap well region due to the viscous 
effects which are dominant In this region. Our 
Inviscid model also Is unable to predict the sep- 
aration that occurs on the upper surface of the 
flap as Indicated by the experimental data. 

Our second computational example deals with 
Iced airfoils. The analysis of the aerodynamic 
performance of Iced airfoils has been of great 
Interest to aircraft designers. In order to find 
ways to prevent Ice formation on wings, one needs 
to accurately predict the flowfleld about Iced 
airfoils with various forms of Ice shapes.*^ 
Ice accretion on an airfoil produces very Irreg- 
ular and rough surfaces on the leading edge 
region. These shapes normally have concavities 
and convexities which cannot be modeled using a 
single zone structured grid. 

An Iced NACA 0012 airfoil with an e1ght-m1nute 
1ce surface computed using the fortified Lewice 
program of Ref, 14 Is shown In Fig. 5. A compos- 
ite structured and unstructured grid Is generated 
consisting of 198 x 32 0-meshes, part of which are 
replaced by 932 triangular meshes. The unstruc- 
tured mesh region Is extended far enough to cover 
the Irregular ice shape on the leading edge, A 
converged solution obtained using this mesh Is 
presented 1n subsequent figures. 

Figure 6 shows the Isobar solution obtained at 
M = 0.2 and a *4®. The compressibility effects 
at this Mach number are less than 2%, according to 
the Prandtl Glauert rule. The contour lines vary 
smoothly across the zonal boundary Indicating that 
the conservation of fluxes Is well satisfied. 

Figure 7 shows the velocity vectors in the 
leading-edge region. This figure clearly Indi- 
cates that there are multiple stagnation regions 
where streamlines approach the surface from the 
farfleld and are divided In two opposite direc- 
tions on the surface. There Is also a small 
region of reverse flow on the upper surface. 
Similar to the first example, although the equa- 
tions are Inviscid, the addition of dissipation 
terms produces viscous flow like vortldty behind 
the Ice shape on the upper surface. 

-The computed pressure distribution on the 
surface Is compared with the surface panel method 
results In Fig. 8; good agreement Is Indicated. 


7.0 Conclusion 

The objective of the present work Is to develop 
an efficient and reliable two-dimensional zonal 
approach capable of coupling structured and un- 
structured grid zones for complex airfoil shapes. 
This approach requires considerably less memory 
compared to using an entirely unstructured mesh 
capable of handling complex geometries, ___It_ also 
takes advantage of the different desirable fea- 
tures of structured and unstructured meshes. 
Computational examples are presented that demon- 
strate the applicability of this method to the 
analysis of the flowflelds such as those around 
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Iced airfoils, multielement airfoils, and blunt 
and divergent-tralling-edge airfoils. The method 
could also be coupled with a boundary-layer method 
to Include the viscous effects. 
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Fig 1 Multielement airfoil with large flap 
deflection. 



solution for multielement airfoil 
of Fig. 1 using composite grid. M = 0.125, a = 0. 



lufnii «l°rV direction for multielement 
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